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ABSTRACT 

A computationally efficient, reduced-order method is presented for prediction of sonic fatigue 
of structures exhibiting geometrically nonlinear response. A procedure to determine the 
nonlinear modal stiffness using commercial finite element codes allows the coupled nonlinear 
equations of motion in physical degrees of freedom to be transformed to a smaller coupled 
system of equations in modal coordinates. The nonlinear modal system is first solved using a 
computationally light equivalent linearization solution to determine if the structure responds 
to the applied loading in a nonlinear fashion. If so, a higher fidelity numerical simulation in 
modal coordinates is undertaken to more accurately determine the nonlinear response. 
Comparisons of displacement and stress response obtained from the reduced-order analyses 
are made with results obtained from numerical simulation in physical degrees-of-freedom. 
Fatigue life predictions from nonlinear modal and physical simulations are made using the 
rainflow cycle counting method in a linear cumulative damage analysis. Results computed 
for a simple beam structure under a random acoustic loading demonstrate the effectiveness of 
the approach and compare favorably with results obtained from the solution in physical 
degrees-of-freedom. 

INTRODUCTION 

Sonic fatigue is an on-going concern for high-speed aircraft and reusable launch vehicle 
structures because of their exposure to high vibroacoustic environments. Such loadings often 
induce geometrically nonlinear response in thin-walled structures, thereby requiring a 
nonlinear random response analysis. Further, a finite-element-based approach is usually 
required due to the geometrical complexity of these structures. Sonic fatigue prediction by 
solution in physical degrees-of-freedom (DoFs) is not attractive because of the computational 
burden associated with a large system size and long simulation time [1] required to obtain 
accurate fatigue estimates. Therefore, reduction of the system is sought to make the solution 
more tractable. 

Within the context of a finite element formulation, several approaches have been developed to 
transform the geometrically nonlinear system in physical DoFs to a smaller coupled nonlinear 
system in modal coordinates. The method employed in this paper indirectly evaluates the 
nonlinear stiffness [2, 3] to allow implementation using the commercial finite element 
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programs MSC.NASTRAN and ABAQUS. Within those codes, standard solution procedures 
are used whenever possible, including pre-processing, eigenanalysis to extract modal basis 
vectors, nonlinear static solutions to determine the nonlinear stiffness, and stress recovery. 
Once the nonlinear stiffness is determined, the reduced-order coupled system of equations is 
solved for the modal displacements, which can subsequently be transformed back to physical 
coordinates and post-processed to obtain stress time histories. 

Two solutions procedures are offered. One is a fast, equivalent linearization procedure, which 
provides useful infonnation about the nonlinear response in an RMS sense. Hence, it is a 
useful tool for screening the structural response for nonlinear behavior. It was previously 
found that equivalent linearization could also provide useful fatigue life estimates, but only 
when the contribution of membrane stress to the total stress was small [4]. The other, more 
robust procedure, is nonlinear modal simulation and is used to provide accurate displacement 
and stress time histories when a nonlinear response is indicated. Response results obtained 
from the reduced-order analyses are compared with those obtained from numerical simulation 
in physical DoFs. Fatigue estimates are computed from stress time histories using a rain flow 
cycle counting analysis and linear cumulative damage. Estimates made using nonlinear 
modal simulation stress time histories are compared with those obtained using results from the 
simulation in physical DoFs. A simple beam structure is used to illustrate the quality of 
reduced-order analysis results achievable, while keeping the numerical simulation analysis 
reasonable. 


NONLINEAR MODAL REDUCTION 

The equations of motion of the nonlinear system in physical DoFs may be written as 

MX {t) + CX{t) + F NL (X(t)) = F{t ) (1) 

where M and C are the system mass and damping matrices, X is the displacement response 
vector and F is the force excitation vector, respectively. As written, the nonlinear restoring 
force F nl is a vector function, which generally includes the linear, quadratic and cubic 
stiffness terms. 

A set of coupled modal equations with reduced DoFs is first obtained by applying the modal 
coordinate transformation X = Oq to Equation (1), where q is the modal displacement 
response vector. The modal basis matrix <D is typically formed from the eigenvectors 
obtained from Equation (1) using only the linear stiffness. Generally, a small set of L basis 
functions are included giving 

Mq(t) + Cq(t) + F nl (q x (t), q 2 ( t ), ...,q L (t)) = F(t) (2) 

where, for mass-normalized eigenvectors, 

M = O t MO = [i] C = O T CO = [2C r 0)J F nl = O t F nl F = O t F (3) 
and co r are the undamped natural frequencies and C, r are the viscous damping factors. 

Nonlinear modal stiffness evaluation. The direct transformation of Equation (1) to (2) 
requires knowledge of the system matrices and excitation vector. In the context of a 
commercial finite element program, the nonlinear restoring force F NL is generally not known. 

Therefore, an indirect means [2, 3] of evaluating the nonlinear stiffness is applied. The 
nonlinear modal restoring force may be written in the form 
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L L L L L L 

^nl id 1 5 • • •’ Ql ) Z^a+EZ zzz*> r = l,2,...,L (4) 

7=1 7=1 k= j 7=1 k= j l=k 

where d, a, and b are the linear, quadratic nonlinear, and cubic nonlinear modal stiffness 
coefficients, respectively. This fonn of F NL reduces the problem of detennining the 

nonlinear stiffness from one in which a large set of simultaneous nonlinear equations must be 
solved to one involving simple algebraic relations. 

The coefficients are obtained through a series of nonlinear static solutions using prescribed 
static displacement fields X c in physical coordinates. For coefficients having all the same 

lower indices, that is, for [fiC] , [a' y ] , and [Zfi ] (j = 1 , 2 ,..., L) , a set of three displacement 
fields are prescribed [3]. For example, for j = 1 , the three displacement fields 

*cl=+4Wl X c2=-hdl X c3 = +4Wl (5) 

are applied. The modal displacements q and q are two different scalar quantities and are 
specified such that the magnitude of the prescribed physical displacement field X c is 

physically meaningful. The corresponding nonlinear restoring force is evaluated using the 
nonlinear static solution within a commercial finite element program, and is given as 
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from which the coefficients [d[\, [u' J and \b r u , ] are obtained. The remaining coefficients 
[d'j J , [a' y ] and [bd] (j = 2,3...,L) are determined in an analogous manner. In this paper, 

modal displacement values of q = 7 x 10 5 and q = 8.75 x 1 0 5 were used for each occurrence. 
It has been shown that the stiffness coefficients obtained via this approach are not sensitive to 
the particular value of modal displacement specified [2, 5]. 

A similar technique can be employed to detennine stiffness coefficients with two unequal 
lower indices, e.g. [a[ 2 \ , [b' n2 \ , and [b r n2 \ . Coefficients of this type appear only if the number 
of retained eigenvectors is greater than or equal to two (Z > 2). Prescribing the displacement 
fields 

X c\ ~ "HWl 4*2^2 X c2 ~ _ 4Wl _ 4*2^2 X c3 ~ ~*~4Wl _ 4*2^2 C^) 

results in the following equations 
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( 8 ) 


from which [<a' 2 ] , [/?,'', 2 ] , and [/?' 22 ] are determined. All remaining coefficients of the type 
[a ' jk ] , [bj jk \ and [b' kkJ \ for j,k = \,2,...,L are found in this manner. 


Finally, for cases when the number of retained eigenvectors if greater than or equal to three 
(L> 3), coefficients with three unequal lower indices, e.g. [/?' 23 ] , are detennined by 
prescribing the displacement field 

X c =+^q l +^ 1 q 2 +^q i (9) 

The resulting equation 

F nl — <D F nl (+(^<7] + §iq 2 + ^,q 2 ) 

— +^1 J ft + \_d 2 J ft ^*^3 J ft 

~*~[^*hl]*7lft "*■ (^**22 J *72*?2 j^ft2 J */lft + [«Jftft [^*^23 *72*73 (10) 

+ [*,'„] ftftft + ^h 2 22 J q 2 q 2 q 2 [^ 333 ] * 73 * 73*73 + [*m] + [«,,] q 2 q 2 q, 

+ [*Tl3 ] ft ft^3 + [^l ] ftftft +[«„] q 2 q 2 q 2 +hi,] q 2 q 2 q 2 + [*m] 


contains one column of unknown coefficients [/ip 3 ] . All remaining coefficients of type 
[bL ] ( / ^ k ^ I) are found in this manner. 


The number of nonlinear static solutions required depends on the number of modes selected 
(L) and may be expressed as: 


Number of Solutions = 3 
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is the number of order-independent combinations of n modes from the L selected modes. For 
example, if L = 12 , the number of nonlinear static solutions required is 454. The stiffness 
evaluation procedure is performed independent of loading. Therefore, the procedure need 
only be performed once prior to the solution of Equation (2). The particular selection of 
modal basis, however, will be influenced by factors including the structural geometry, 
material properties, and loading distribution. Recent work by the authors investigated basis 
selection for flat [5] and curved structures [6]. 
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Implementation. The nonlinear stiffness evaluation method, and numerical integration and 
equivalent linearization solutions to follow, are embodied in the program RANSTEP, for 
“Reduced order Analysis using a Nonlinear STiffness Evaluation Procedure.” RANSTEP has 
been implemented for use with two commercial finite element programs; MSC.NASTRAN 
version 2005 and ABAQUS version 6.5. 

To obtain the modal basis matrix <D for the coordinate transformation X = Q)q , the 
NASTRAN implementation uses a direct matrix abstraction program (DMAP) to write mass- 
normalized eigenvectors from the standard normal modes analysis (solution 103). The 
ABAQUS implementation uses a PYTHON script to access mass-normalized eigenvectors 
from the output database created using the *FREQUENCY analysis in ABAQUS Standard. 

For both implementations, an external FORTRAN program forms the set of prescribed 
displacement fields, as per Equations (5), (7) and (9). Since the displacement fields are 
fonned outside of the finite element program, other basis functions may be utilized, e.g. 
experimentally determined, however these are not considered in the present work. A static 
nonlinear analysis follows to determine the corresponding nonlinear restoring forces. In the 
NASTRAN implementation, a DMAP code is used to modify the standard nonlinear static 
analysis (solution 106) to allow specification of the displacement fields and output of the 
nonlinear restoring forces. In the ABAQUS implementation, a PYTHON script extracts the 
forces from the output database created using the * STATIC analysis (with geometric 
nonlinearity enabled) in ABAQUS Standard. The nonlinear stiffness coefficients are 
determined from Equations (6), (8) and (10) in an external FORTRAN code using the 
nonlinear forces from the above analyses. 

Nonlinear stiffness validation. Validation of the nonlinear modal stiffness fonnulation and 
implementation was performed using an aluminum beam under a uniformly distributed 
pressure loading. The beam measured 18-in. x 1-in. x 0.09-in (/ x w x h ) and was fully 
clamped at both ends. The following material properties were used: 

E = 10.6 x \0 6 psi, G = 4.0 x \0 6 psi, p = 2.588 x 10^ 

Displacements and stresses obtained from nonlinear modal solutions using different modal 
bases were compared with those obtained directly from a full nonlinear analysis in physical 
DoFs. A static loading was selected to eliminate any dynamic effects associated with the 
modal mass and damping. Hence, only the modal stiffness coefficients dictate the accuracy of 
the solution for a particular basis. While it is possible to set the modal accelerations and 
velocities to zero in Equation (2) and solve for the modal displacements via Newton-Raphson 
(as per [2]), the more direct approach utilized here applies the load quasi-statically and 
numerically integrates Equation (2) to obtain the modal displacements. The modal numerical 
integration procedure used is discussed in the next section. The modal displacements are 
transformed to physical coordinates and stresses are obtained through post-processing, in a 
procedure also discussed in the next section. 

The quasi-static pressure loading was applied to the beam with a loading time history ramped 
from zero magnitude at time zero to its maximum level at 1.6384s. The response was 
obtained after an additional 0.5s at the maximum loading level to allow for the decay of any 
transient behavior. 

For the reduced order analysis, the beam was modeled in MSC.NASTRAN using 144 
CBEAM elements. The analysis in physical DoFs was perfonned using the ABAQUS 
explicit solution. The double precision explicit integration scheme with an adaptive time 
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integration step (referred to as “element by element” in ABAQUS) was utilized for the quasi- 
static and subsequent dynamic analyses. The ABAQUS model consisted of 144 B21 
elements. Differences between the ABAQUS and MSC.NASTRAN element fonnulations 
were previously found to be negligible [5]. 

It was previously shown that a basis consisting of a single mode could exactly reconstruct the 
displacement when the applied loading distribution was identical to the mode shape [2], For 
other distributions, bases consisting of a number of modes are required to obtain an accurate 
solution. For the present analysis, two modal bases were considered. The first basis consisted 
of a combination of uncoupled low-frequency bending modes and high-frequency membrane 
modes. The inclusion of both bending and membrane modes in the basis captures the effect 
of nonlinear bending-membrane coupling. In particular, the basis considered was comprised 
of the first six symmetric bending modes and the first six anti-symmetric membrane modes. 
A second modal basis consisting of only the first six bending modes was also considered. 
This basis is consistent with subsequent linear and equivalent linear analyses, for which there 
is no nonlinear bending-membrane coupling. As such, the basis is not capable of predicting 
the membrane displacement response under a transverse loading. 

Shown in Figure 1 is the transverse displacement response for the two modal bases. The 
displacement results obtained using both bases are comparable to each other and compare 
well with the physical DoF solution As expected, only the basis consisting of bending and 
membrane modes compares well with the physical DoF solution for the membrane 
displacement, as shown in Figure 2. The stress response is shown in Figure 3 and Figure 4 for 
the bending and membrane components, respectively. The fact that there are several points 
along the beam where membrane stress results from the bending only basis agrees with the 
physical DoFs solution highlights the need to consider a number of spatial locations in 
assessing the accuracy of the modal basis. 




Figure 1: Quasi-static transverse Figure 2: Quasi-static membrane 

displacement along beam under 0.0576 psi displacement along beam under 0.0576 psi 

pressure load. pressure load. 
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Figure 3: Quasi-static bending stress along Figure 4: Quasi-static membrane stress 
beam under 0. 05 76 psi pressure load. along beam under 0. 05 76 psi pressure load. 


A compact measure of the displacement and stress error along the span of the beam, relative 
to the physical DoF solution, may be found from 


% Error = 



(Physical i - Reduced Order 'f 
Physical ] 


xlOO 


(13) 


where A is a number of points along the beam span used for the estimation. Figure 5 and 
Figure 6 present the displacement and stress errors, respectively, computed using results 
obtained at the seven locations (N= 7) shown in Figure 1 - Figure 4. It is seen that the bending 
and membrane basis provides excellent agreement through the entire range of excitation 
levels studied, and validates both the nonlinear modal stiffness formulation and the 
implementation. The large error associated with the bending only basis is indicative of the 
quality of results that may be achieved in the subsequent linear and equivalent linear analyses. 
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Figure 5: Quasi-static displacement error 
integrated along the beam semi-span. 


Figure 6: Quasi-static stress error 
integrated along the beam semi-span. 
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SOLUTION PROCEDURES 


Two solution procedures of differing fidelity and computational burden are considered. The 
equivalent linearization solution is fast, but accurate only in an RMS sense. The nonlinear 
modal simulation is accurate in terms of its ability to capture nonlinear spring stiffening and 
softening, and broadening of the power spectral density. It comes at greater computational 
expense relative to equivalent linearization, but less relative to simulation in physical DoFs. 

Equivalent linearization. An approximate solution to Equation (2) can be achieved by 
formation of an equivalent linear system: 

Mq(t) + Cq(t) + K eq q(t) = F(t) (14) 

where K eq is the modal equivalent linear stiffness matrix given by 

K eq =O r K eq O. (15) 


Using the force-error minimization method of equivalent linearization, the modal equivalent 
linear stiffness matrix is determined in an iterative approach, as discussed in [7]. For the first 
iteration, K eq is simply the linear modal stiffness coefficient matrix. For subsequent 
iterations (m> 1), the nonlinear stiffness is written as 


d \ + X a \ E [ q j ] + X X h )n E \ q flk ] 

7 j I k j 

X a n E [ q j ] + X X b jki E [ q j q k ] 

j 7=1 k=j 


X a \ E [ q j ] + X X b )kL E [ q j q k ] 

j 7=1 A j 

d L + X a 'ji ^ M + X X b k E [ q j q k ] 

j 7=1 k=j 


Once K eq is known, Equation (14) can be solved in the time domain or frequency domain, 

with subsequent modal results transformed back to physical coordinates. In this paper, the 
equivalent linear system is numerically integrated so that response time histories may be 
obtained. 

Nonlinear modal simulation. For the nonlinear modal simulation procedure, the coupled 
modal nonlinear equations of motion in Equation (2) are numerically integrated using a 
fourth-order Runge-Kutta method. The resulting modal displacement time histories are 
transformed back to physical coordinates using the inverse modal transformation. 

The equivalent linearization and numerical simulation solutions have been implemented in 
FORTRAN programs in RANSTEP, and operate independent of the particular finite element 
program used to generate the modal nonlinear stiffness coefficients. 

Stress recovery. Element stress recovery is perfonned by post-processing the nodal physical 
DoFs directly within the finite element program. A stress time history is obtained by 
specifying the physical DoFs for the element(s) of interest at each simulation output time step. 
As a result, the stresses are identical to those that would have been obtained via a finite 
element analysis in physical DoFs. 

For the equivalent linear analysis, stress time histories are recovered using a linear static 
analysis procedure with specified nodal displacements obtained from each integration output 
time step. In MSC.NASTRAN, the linear static solution 101 with multiple subcases is used. 
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In ABAQUS Standard, the ^STATIC solution is used. For either finite element program, the 
stress recovery mesh need only encompass the element or elements of interest, and therefore 
the process generally has a light computational burden. A nuance of MSC.NASTRAN 
requires the additional attachment of a dummy element and load to the stress recovery 
element, in order to calculate the resulting element nodal forces and recover the element 
stress. 

For the nonlinear modal simulation analysis, stress time histories are instead recovered using 
a nonlinear static analysis procedure. In MSC.NASTRAN, the nonlinear static solution 106 
with multiple subcases is used. In ABAQUS, either the Standard analysis (using the 
*STATIC solution with geometric nonlinearity enabled) or Explicit analysis is used. The 
ABAQUS Explicit analysis has been found to be significantly faster than the Standard 
analysis, but requires that the initial conditions additionally be specified in the first analysis 
step. This is because the Explicit analysis is still being run as if it were a dynamic analysis 
and not a series of unrelated static cases. 

RANDOM RESPONSE 

A clamped-clamped beam under a uniformly distributed random pressure loading with a 
bandwidth of 1500 Hz was chosen to demonstrate the quality of results obtainable via a linear 
modal solution, equivalent linear solution, nonlinear modal simulation, and numerical 
simulation in physical DoFs. This structure was selected as it exhibits nonlinear bending- 
membrane coupling and is sufficiently small in size to perfonn the numerical simulation in 
physical DoFs within a reasonable period of time. The analysis of a more practical structure 
is the subject of reference [8]. 

The linear modal, equivalent linear and nonlinear modal simulation analyses use the 
MSC.NASTRAN beam model and properties previously discussed. The linear modal and 
equivalent linear analyses used the lowest six symmetric bending modes in the modal 
transformation, while the nonlinear modal simulation used both the lowest six symmetric 
bending modes and the lowest six anti- symmetric membrane modes. The numerical 
simulation in physical DoFs was performed using the ABAQUS explicit solution and the 
beam model previously discussed. Damping was chosen to be sufficiently high so that good 
comparisons could be made at the peaks of the power spectral density (PSD) function. A 
level of mass proportional damping was specified corresponding to 2.0% critical damping for 
the first symmetric bending mode. 

Load generation and ensemble averaging. The same loading time history and ensemble 
averaging was used for all analyses. The loading time histories were generated by summing 
equal amplitude sine waves, each with random phase, within the specified bandwidth using a 
discrete inverse Fourier transform. This procedure was identical to that used in previous work 
[7] by the authors, so further details are omitted for brevity. The loading produced by this 
method has a Gaussian distribution. A sharp roll-off of the input spectrum practically 
eliminates excitation of the structure outside the frequency range of interest. 

A range of load levels was considered to span the response regime from nearly linear to 
highly nonlinear. For each load level, ten ensembles of displacement and stress response 
were generated lasting 2.138% each. The first 0.5s of each response record was discarded to 
remove the initial transient response [7], resulting in response histories of 1.6384s in duration. 
For each simulation, the displacement and stress response was stored at every 50 //s , giving 
time records of 32,768 points. A 32,768-point FFT was subsequently used to compute the 
PSD function. 
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Results. Displacement and stress results are presented at approximately the %-span location 
to highlight the contribution of the membrane response. Displacement results were obtained 
at the node located at 4. 5 -in. from the clamped end. Element stress results were obtained at 
an adjacent element, with a stress recovery point at 4.4375-in. from the clamped end. 

The root-mean-square (RMS) transverse displacement is shown in Figure 7 for a range of 
loadings. The equivalent linear analysis indicates a nonlinear response. Thus, there is a need 
to run the more accurate nonlinear modal simulation to obtain stress time histories needed for 
the fatigue estimate. The best comparison with numerical simulation in physical DoFs is 
achieved using the nonlinear modal simulation, which slightly over-predicts the response. 
Good agreement is also achieved, however, with the equivalent linear analysis. Note that the 
inclusion of six additional anti- symmetric membrane modes (not shown) in either the linear 
modal or equivalent linear analyses does not significantly change the quality of the 
comparison. This is due to the lack of linear bending-membrane coupling for this problem. 
The RMS membrane displacement is shown in Figure 8. Compared to the physical DoF 
results, the level of agreement achieved using the nonlinear modal simulation is very good. 
Since the membrane displacement only manifests itself via nonlinear bending-membrane 
coupling for this problem, the linear modal and equivalent linear analysis are incapable of 
generating any membrane displacement component for a basis consisting of only bending 
modes. 



Figure 7: RMS transverse displacement as 
a function of random pressure load. 
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Figure 8: RMS membrane displacement as 
a function of random pressure load. 


The transverse displacement PSD for the highest RMS pressure of 0.4608 psi (164 dB) is 
shown in Figure 9. For this level, the response is highly nonlinear and exhibits nonlinear 
spring hardening, seen by a shift in the peak response to higher frequencies relative to the 
linear analysis and peak broadening. Very good agreement between the nonlinear modal and 
physical DoFs solution is seen across the frequency range. The equivalent linear analysis, 
which was seen to agree fairly well in an RMS sense, is significantly different from the 
physical DoFs solution, highlighting the need for a more accurate response calculation. The 
magnitude of the peak response is different, the shifting of peaks to higher frequency is 
exaggerated, and it lacks the peak broadening characteristic. The linear modal response, as 
expected, captures none of the essential characteristics. The corresponding membrane 
displacement PSD, shown in Figure 10, indicates that the nonlinear modal simulation captures 
the essential features of period doubling, with greater differences observed at the higher 
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frequencies. Again, the linear modal and equivalent linear analyses are incapable of 
predicting the membrane displacement. 




Figure 9: Transverse displacement PSD for Figure 10: Membrane displacement PSD 

a 164 dB random pressure loading. for a 164 dB random pressure loading. 


Figure 1 1 shows the normal stress PSD on the beam surface corresponding to the highest 
loading of 164 dB. The nonlinear modal simulation agrees very well with the solution in 
physical DoFs. Note the significant zero frequency component in both the nonlinear modal 
simulation and physical DoFs analyses. This is due to the membrane stress component, which 
oscillates between zero and some positive value. For the physical and nonlinear modal 
simulations, the membrane stress skews the probability distribution function (PDF) to the 
positive stress range, as shown in Figure 12. Flere the PDF has been normalized by the 
standard deviation. The equivalent linear and linear modal stress PSDs bear little 
resemblance to the solution in physical DoFs. Further, these analyses are unable to predict the 
membrane stress. Consequently, the zero frequency component is small for the equivalent 
linear and linear modal solutions and their PDFs are nearly Gaussian, with zero skewness. 



random pressure loading. 



Figure 12: Normalized normal stress PDF 
for a 164 dB random pressure loading. 


To recap, the equivalent linear analysis was seen to be a good indicator of nonlinear response 
through observation of the RMS transverse displacement. As expected, it was not sufficient 
to adequately compute the membrane behavior, the displacement PSD, and stress PSD. 
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Hence, a more accurate means of computing these quantities is required. The nonlinear 
modal simulation was able to meet that need, provided a sufficient modal basis is used. 
Therefore, in the following fatigue analysis, results from the nonlinear modal simulation are 
considered as the only alternative to the results from the simulation in physical DoFs. 

FATIGUE ANALYSIS 

For a constant amplitude, zero-mean stress response, the stress-life (S-N) curve relates the 
cycles to failure, N f , to the fully reversed alternating stress, <r ar , via the power law: 


^ AN* 


(16) 


where A = cr'. 2 ’ , and cr' f and b are material properties determined from tests perfonned under 
zero-mean stress. The fatigue life can then be determined from 


N f = 


f (j \/b 


v A J 


(17) 


For conditions having a non-zero mean stress, a number of models are available, see [9], In 
this paper, we consider the Walker model, which expresses the alternating stress as: 
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and the number of cycles to failure as: 
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(19) 


In the above, a' fw , b w y are material properties determined from tests performed under a 
range of stress ratios R. For fully reversed loadings, R = -1 and <J ar = cr max . 


The Palmgren-Miner linear cumulative damage rule [10] is typically used for variable 
amplitude loading and assumes that the damage, D, caused by stress cycles at one stress level 
can be calculated and added to damage caused by stress cycles in another stress level, or 


"=Z 
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( 20 ) 


where IV, are the number of cycles at stress level cr. , and (/V^j are the number of cycles to 
failure at stress level cr . 


Random fatigue. For random response, it is convenient to recast Equation (20) in the form 


e[d] = e[p]tY^ 


p{& ml n^max) A °'min A(J lr 
Nf (^min ’ ^max ) 


( 21 ) 


where E[P ] are the number of peaks per second, T is the lifetime, />(cr min ,cr max ) is the joint 
PDF of stress minimum and maximum values, and Act . and A cr are the PDF bin widths. 

J mm max 
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The dependency of N f on cr min and cr max is detennined by the particular model used, e.g. 
zero-mean. Walker, etc. 


Rainflow analysis. It was shown in [11] that p(cr min , cr mx ) can be estimated from the 
rainflow matrix (RFM) 

J oo oo oo oo 

t^EE RFM a msx ) \ Tnin ’ '-’max ) ^ Tniri ^ ^nm (22) 

1 ' Rf -00 —00 —00 — OO 

where N RF are the number of rainflow cycles. The RFM is a two-dimensional histogram of 
rainflow ranges. In this paper, the RFM was obtained using the WAFO Matlab toolbox for 
analysis of random waves and loads [12]. For a finite record of duration 7] nax , 

E[P] = N rf / T max . Substitution of Equation (22) into (2 1 ) then yields 


£[£>] 


T y. y. flFM (tT ||lin , tT |nax ) 
17 LuLu A T ^ \ 


' T ’ 00 00 

'V 'V RFD(a ,cr ) 

ry-7 / ^ V min 9 max / 

max _co _c0 


where RFD(a min ,a max ) is the two-dimensional rainflow damage matrix. Failure occurs when 
E[D\ = 1 , giving the estimated fatigue life T as 


T = 


/ 00 00 

EE 


RFM (a ,o" ) 

V mm 9 max / 

(^min ’ '-’"max ) 


(23) 


Note that since damage is calculated for each bin of the RFM, the quantities in Equation (18) 
vary with each cycle. 


Fatigue results. Fatigue life estimates were made at two locations on the beam previously 
considered. A position near the clamped end (0.0625-in. from the clamped end) was selected 
because it exhibits the highest stress and has the shortest fatigue life. A second position near 
the 14-span (4. 43 75 -in. from the clamped end) was selected to show the effect of a significant 
membrane component on the fatigue life prediction. At each loading level, ten normal stress 
history records were joined together to form an overall record length of T max = 16.3845. S-N 
properties for 2024-T3 aluminum were used [9], with cr' f = 1602 and b = -0.154 for the zero 

mean case, and cr' /u =1772, b w = -0.163, and y = 0.460 for the Walker model, when the 

stress is expressed in units of MPa. For the Walker fatigue estimate, damage was not 
accumulated for cycles with R<- 2 as the S-N data is not considered valid in this range. 


It is useful to first investigate the rainflow matrix. The RFM for the clamped end at the 164 
dB loading is shown in Figure 13. The RFM in the upper left domain was generated from 
the stress time history obtained from the simulation in physical DoFs. Superimposed in the 
lower right domain is the RFM obtained from the nonlinear modal simulation. If the two 
were identical, they would form the reverse mirror image of each other. Recall that the 
membrane stress component accounts for the zero frequency component and hence for the 
skewness in the normal stress response. Even though the normal stress is dominated by the 
bending component at this location, the membrane stress still skews the RFC matrices, as 
evidenced in data from both for the physical and nonlinear modal simulations. The solid 
diagonal line in Figure 13 indicates the fully-reversed (R - -1) rainflow cycles. 
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Fatigue life estimates at the clamped end are presented in Table 1. Estimates were made 
using the zero-mean and Walker models. The estimates are expressed in ranges and were 
obtained by scaling the stress time histories up and down by 10%. This was done in an effort 
to reflect the sensitivity of the life prediction to variability in the stress prediction. Also 
shown are the stress standard deviation and mean. Previous work indicated that zero-mean 
models produce non-conservative fatigue estimates relative to non-zero mean models for 
tensile mean stress conditions [11, 13]. Estimates made using data from both physical and 
nonlinear modal simulations confirm this finding. It is interesting to note, however, that as 
the mean stress component increases from roughly 5% of the standard deviation at the 146 dB 
level, to roughly 14% at the 164 dB level, the difference between zero-mean and non-zero 
mean life estimates is reduced from about a factor of 3 (153 vs 53 years, for example), to less 
than a factor of 2 (3.05 vs 1.68 hours). In general, life estimates made using nonlinear modal 
simulation stress time histories compare favorably, though somewhat non-conservatively, 
with those made using physical simulation stress time histories. To understand why, it is 
helpful to look at the rainflow damage matrix. Figure 14 shows the Walker rainflow damage 
matrices corresponding to Figure 13. The RFD in the upper left domain was generated from 
the physical simulation RFM, while the RFD superimposed in the lower right domain was 
generated from the nonlinear modal simulation RFM. More damage occurs at the high stress 
amplitudes for the physical DoFs solution than at the high stress amplitudes for the nonlinear 
modal simulation. Consequently, life estimates made from physical simulation data are lower 
than those obtained using nonlinear modal data. 



Figure 13: Rainflow matrices at 
clamped end for 164 dB loading. 



Min Stress Amplitude (ksi) 


Figure 14: Walker rainflow damage 
matrices at clamped end for 164 dB loading. 


Table 2 presents fatigue life estimates at the Vi-span location, where the membrane stress is 
more significant. While this location will not dictate the failure of the beam structure, it is 
nevertheless interesting to note that fatigue life estimates from the nonlinear modal simulation 
agree very well with those obtained from the physical DoFs simulation. 


In summary, both zero-mean and Walker fatigue life estimates obtained using nonlinear 
modal simulation data compare quite favorably with estimates made using physical simulation 
data, irrespective of the location, because the stress response from which the estimate is 
derived is in close agreement. 
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Table 1: Stress standard deviation, mean, and fatigue life estimates 
at the clamped end for a range of loadings. 



Stress Std. Dev. / Mean 
(psi) 

Zero Mean 
Fatigue Fife Estimate 

Walker 

Fatigue Fife Estimate 

OASPL 

(dB) 

Physical 

DoFs 

Nonlinear 

Modal 

Physical 

DoFs 

Nonlinear 

Modal 

Physical 

DoFs 

Nonlinear 

Modal 

146 

1379/69 

1332/68 

153-563 yr 

215-791 yr 

53.0-181 yr 

73.8-253 yr 

152 

2685 /210 

2466/ 194 

1.65-6.08 yr 

3.07-11.3 yr 

0.63-2.17 yr 

1.24M.28 yr 

164 

8751 / 1211 

8235 / 1170 

3.05-11.2 hr 

6.17-22.7 hr 

1.68-5.74 hr 

3.49-12.0 hr 


Table 2: Stress standard deviation, mean, and fatigue life estimates 
at the %-span for a range of loadings. 



Stress Std. Dev. / Mean 
(psi) 

Zero Mean 
Fatigue Fife Estimate 

Walker 

Fatigue Fife Estimate 

OASPF 

(dB) 

Physical 

DoFs 

Nonlinear 

Modal 

Physical 

DoFs 

Nonlinear 

Modal 

Physical 

DoFs 

Nonlinear 

Modal 

152 

868/215 

815/198 

1823-6711 yr 

2796-10289 yr 

387-1324 yr 

554-1897 yr 

164 

2986 / 1224 

2951 / 1269 

0.61-2.24 yr 

0.72-2.64 yr 

0.13-0.47 yr 

0.15-0.50 yr 


CONCLUDING REMARKS 

Two reduced-order analysis methods, equivalent linearization and nonlinear modal 
simulation, have been presented. These methods have been implemented in the program 
RANSTEP, which utilize the finite element programs MSC.NASTRAN and ABAQUS to 
evaluate the nonlinear modal stiffness and perform stress recovery. The nonlinear modal 
stiffness formulation and implementation has been shown to recover the solution in physical 
degrees-of-freedom when a suitable modal basis is used. 

The two analysis methods were shown to offer complementary capabilities for the analysis of 
high cycle fatigue. The equivalent linearization analysis was shown to be an effective 
screening tool to detennine if a higher fidelity nonlinear analysis is required. If no significant 
nonlinearity is indicated, then fatigue estimates based on the customary linear analysis should 
be adequate. Similarly, if significant nonlinearity is indicated, a nonlinear modal simulation 
may be used to efficiently compute the fatigue life. Since the accuracy of either reduced- 
order method is determined largely by the modal basis used, confidence in the results obtained 
must be established through comparison with either experimental data or other representative 
analyses. In the present work, displacement response, stress response, and fatigue life 
predictions made using nonlinear modal simulation results compared favorably with those 
obtained from the solution in physical DoFs. 
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